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We provide a phenomenological theory for topological transitions in restructuring networks. In this statistical 
mechanical approach energy is assigned to the different network topologies and temperature is used as a quantity 
referring to the level of noise during the rewiring of the edges. The associated microscopic dynamics satisfies the 
detailed balance condition and is equivalent to a lattice gas model on the edge-dual graph of a fully connected 
network. In our studies - based on an exact enumeration method, Monte-Carlo simulations, and theoretical 
considerations - we find a rich variety of topological phase transitions when the temperature is varied. These 
transitions signal singular changes in the essential features of the global structure of the network. Depending 
on the energy function chosen, the observed transitions can be best monitored using the order parameters $ s = 
Smax/Af, i.e., the size of the largest connected component divided by the number of edges, or <E>fc = fc max /M, 
the largest degree in the network divided by the number of edges. If, for example the energy is chosen to be 
E = — Smax, the observed transition is analogous to the percolation phase transition of random graphs. For 
this choice of the energy, the phase-diagram in the [{k) , T] plane is constructed. Single vertex energies of the 
form E — £\ f(ki), where ki is the degree of vertex z, are also studied. Depending on the form of f(ki), first 
order and continuous phase transitions can be observed. In case of /(fcj) = — (ki + a) ln(fci), the transition is 
continuous, and at the critical temperature scale-free graphs can be recovered. Finally, by abruptly decreasing 
the temperature, non-equilibrium processes (e.g., nucleation, growth of particular topological phases) can also 
be interpreted by the present approach. 

PACS numbers: 9.75.Hc, 05.70.Fh, 64.60.Cn, 87.23.Ge 



I. INTRODUCTION 

In recent years, the analysis of the network structure of in- 
teractions has become a popular and fruitful method used in 
the study of complex systems. Whenever many similar objects 
in mutual interactions are encountered, these objects can be 
represented as nodes and the interactions as links between the 
nodes, defining a network. The world-wide-web, the science 
citation index, and biochemical reaction pathways in living 
cells are all good examples of complex systems widely mod- 
eled with networks, and the set of further phenomena where 
the network approach can be used is even more diverse. In 
most cases, the overall structure of networks reflect the char- 
acteristic properties of the original systems, and enable one to 
sort seemingly very different systems into a few major classes 
of stochastic graphs 1 1, 2\. These developments have greatly 
advanced the potential to interpret the fundamental common 
features of such diverse systems as social groups, technolog- 
ical, biological and other networks. The effects of both the 
restructuring |3] and the growth |4] of the associated graphs 
have been considered, leading to a number of exciting dis- 
coveries about the laws concerning their diameter, clustering 
and degree distribution. Real networks typically exhibit both 
aspects (growth and rearrangement) one of which is usually 
dominating the dynamics. Here we concentrate on the evo- 
lution of graphs due to restructuring, but shall briefly discuss 
the growth regime as well. 

Various interesting effects observed in networks can be in- 
terpreted using analogies with well understood phenomena 
studied in statistical physics. As a classical example, we men- 
tion the percolation phase transition in the Erdos-Renyi (ER) 



random graph model |5l|CJ], which occurs by varying the aver- 
age degree, (k), of the vertices around (k) = 1. For (k) < 1 
the graph falls apart into small pieces, on the other hand for 
(k) > 1 a giant connected component emerges (in addition to 
the finite components). Another subtle example is the map- 
ping of a growing network model onto an equilibrium Bose 
gas [7]. For the latter model, under certain conditions a single 
node is allowed to collect a finite fraction of all edges, corre- 
sponding to a highly populated ground level and sparsely pop- 
ulated higher energies seen in Bose-Einstein condensation. 

When connecting the graph theoretical aspects of networks 
to statistical physics, one can step further from the analogies 
by directly defining statistical ensembles for graphs. The use 
of a statistical mechanical formalism for the changes in graphs 
being in an equilibrium-like state is expected to provide a sig- 
nificantly deeper insight into the processes taking place in sys- 
tems being in a saturated state and, as such, dominated by the 
fluctuating rearrangements of links between their units. 

As an example, let us take a given number of units inter- 
acting in a "noisy" environment. These units can be people, 
firms, genes, etc. The probability for establishing a new or 
ceasing an existing interaction between two units depends on 
both the noise and the advantage gained (or lost) when adopt- 
ing the new configuration. In this picture, a global transition 
in the connectivity properties can occur as a function of the 
level of noise. For instance, if the conditions are such that 
the interactions between the partners become more "conser- 
vative" (safer choices are more highly valued), then - as we 
show later - a transition from a less ordered to a more ordered 
network configuration can take place. In particular, it has been 
argued |8] that depending on the level of certain types of un- 



2 



certainties (expected fluctuations) business networks reorga- 
nize from a star-like topology to a system of more cohesive, 
highly clustered ties. 

There are several possible ways to define the statistical en- 
semble of networks. In I9t llOII . the members of the ensem- 
ble were identified by the Feynman diagrams of a field theory 
in zero dimensions (called "minifield"), and the weights of 
the graphs were given by the corresponding amplitudes cal- 
culated using the standard Feynman rules. The ensemble ob- 
tained this way was characterized by the fractal and spectral 
dimensions, and the dependence of the topology of the graphs 
on these two parameters was discussed. The authors argued 
that in the parameter plane of two parameters related to the 
fractal and spectral dimension, the region of generic graphs 
and the region of crumpled graphs are separated by a line; 
and on this separating line scale-free networks appear. An 
alternative definition for the partition function was proposed 
by Berg and Lassig in full , resulting in a simpler formalism, 
analogous to the statistical mechanics of classical Hamilto- 
nian systems. They introduced a Hamiltonian for networks, 
and also a parameter (3 playing the role of inverse temperature. 
The weights of different graphs in the partition function were 
obtained from these two quantities as in classical statistical 
mechanics. These studies showed that Hamiltonians beyond 
the single vertex form (where terms depending on connectivi- 
ties between the vertices also appear) lead to correlations be- 
tween the vertices for large (3. A similar model leading to 
interesting results was presented in where the Hamil- 
tonian depended on the ratios of the degrees of neighboring 
vertices, and the dynamics favored disassortative mixing and 
high clustering. The system organized itself into three phases 
depending on one parameter: the exponential, scale-free and 
hub-leaves states were produced, respectively. However, the 
non-uniform selection of links at the rewiring in this model 
makes it impossible to satisfy the detailed balance condition. 

In this paper we analyze the reorganization of networks 
from the point of view of topological phase transitions, i.e., 
transitions in the graph structure as a function of tempera- 
ture, the quantity representing the level of noise during the 
restructuring process of the network. For clarity we note that 
our studies concern a class of phenomena that are clearly 
different from the phase transitions investigated by applying 
models of statistical mechanics originally defined on regu- 
lar lattices to an underlying (static) random network struc- 
ture Il 4ill 3|. or the phase transitions observed in growing 
networks |7, 16], or quasi-static networks |17]. Topological 
phase transitions are accompanied by singularities in the ther- 
modynamic functions derived from the partition function of 
the statistical graph ensemble and can be characterized by a 
drastic change in an appropriate order parameter. Our statis- 
tical ensemble (similarly to the one presented in Ref. flllo . is 
defined by introducing an energy that accounts for the advan- 
tage or loss during the rearrangement. In our description, this 
energy may depend on either global properties of the network, 
or single vertex degrees as well. 

The use of Hamiltonian formalism also provides a general 



frame for the optimization of network structure (for examples 
of network optimization problems see I18lll9l0 . To find the 
optimal configuration for a given task, the system has to be 
cooled, using an appropriate energy function. 

This article is a direct extension of our previous work 1 20] ; 
covering more details, results, and new approaches. The pa- 
per is organized as follows: in Sec. II we define the canonical 
ensemble of the networks together with the partition function 
and other essential thermodynamic quantities. In Sec. Ill we 
discuss the numerical methods used to study the phase tran- 
sitions. In Sec. IV we present the phase transitions obtained 
for energies that depend on global properties, and Sec. V is 
devoted to two interesting cases of single vertex energies. In 
Sec. VI we discuss briefly the grand canonical ensemble of 
networks and we conclude in Sec. VII. 



II. STATISTICAL MECHANICS OF NETWORKS 

We shall consider a set, {g a }, of undirected graphs, con- 
taining N nodes and M links. Each graph g a can be repre- 
sented by the adjacency matrix Afj, where A!L — 1 if vertices 
i and j are connected and it is zero otherwise. In a heat bath 
at temperature T, the canonical ensemble of these graphs (in 
analogy with that proposed in HHl "). can be defined by the par- 
tition function 

Z(T) = £ c- E °/ T , (1) 

{9a} 

where E a is the energy assigned to the different configura- 
tions. 

The restructuring processes of the network can be inter- 
preted via the following physical picture: The basic event 
of rearrangement is the reallocation of a randomly selected 
edge (link) to a new position either by "diffusion" (keeping 
one end of the edge fixed and connecting the other one with 
a new node) or by removing the given edge and connecting 
two randomly selected nodes. Then, the energy difference 
AE a b = Et — E a between the original g a and the new gi, 
configurations is calculated and the reallocation is carried out 
following the Metropolis algorithm 1 2 1 ] . If the energy of the 
new graph is lower than that of the original one, the reallo- 
cation is accepted; if the new energy is higher, the realloca- 
tion is accepted only with probability e - AE ab/T _ wa y^ m 
the T — > oo limit the dynamics converges to a totally random 
rewiring process, and thus, the classical ER random graphs are 
recovered. On the other hand, at low temperatures the topolo- 
gies with lowest energy occur with enhanced probability. The 
resulting dynamics, by construction, satisfies the detailed bal- 
ance condition I2H1 . 

This network rearrangement is formally equivalent to a 
Kawasaki type lattice gas dynamics with conserved number 
of particles moving on a special lattice,which is the edge-dual 
graph of the fully connected network |6, 22]. The sites of this 
lattice are the possible N{N — l)/2 connections between the 



3 






FIG. 1: Two simple examples of graphs (left hand side) and the 
corresponding edge-dual graphs (right hand side). The full spheres 
in the edge-dual graphs represent occupied sites, corresponding to 
existing bonds in the original graph (drawn with solid lines), and the 
hollow spheres are the empty sites, corresponding to absent bonds 
(represented by dashed lines) in the original network. The rewiring 
of an edge in the original graphs is equivalent to the displacement of 
the corresponding particle on the edge-dual graph. 

vertices, and the particles wandering on the sites are the M 
edges, as shown in Fig. [2 

The partition function Q contains many terms correspond- 
ing to topologically equivalent graphs: these graphs can be 
simply transformed into one another by an adequate permu- 
tation of the indexing. Since we consider energies E a that 
depend only on the topology t a , it is natural to rewrite the par- 
tition function in a form where the summation runs through all 
possible topologies: 



Z(T) = ^A/;e- B =/ T . 



(2) 



Here we introduced Af a to count the number of configurations 
belonging to topology t a . Expression (0 can be rewritten as 



Z(T) = J2 e- E ^ T+ln ^ = e"^ /T , 



{*<*} 
S a = \n{M a ), 



(3) 

(4) 
(5) 



where F a is the free energy and S a is the entropy of the topol- 
ogy t a . 

We are interested in the possible singularities in the thermo- 
dynamic functions derived from the partition function above, 
since they, if there are any, correspond to phase transitions 
in the topology of the associated networks. These transitions 
can be best monitored by introducing a suitable order param- 
eter. As we are primarily interested in the transitions between 
dispersed and compact states, a natural choice can be either 



$ = $ s = s max /M, the number of edges of the largest con- 
nected component of the graph s max normalized by the total 
number of edges M, or $ = $fe = k max /M, the highest de- 
gree in the graph /c max divided by M. We also introduce the 
corresponding conditional free energy F(<!>, T) via 



e- F W r >/ r = Z(*,T)= Y, e ~ £ " /T ' 



(6) 



where {</ a }$ is a subset of {g a }, consisting of all the graphs 
with order parameter $. A phase transition, where a rapid 
change occurs in the order parameter from $ = towards 
higher values, is also accompanied by a shift of the minimum 
of the conditional free energy F(<£>, T). A sudden change in 
the position of the global minimum signals a discontinuous 
(first order) phase transition, whereas a gradual shift indicates 
either a crossover or a continuous phase transition. 

In the next section we briefly discuss the numerical methods 
used to study topological phase transitions. 



III. NUMERICAL METHODS 

///. A ) Exact enumeration method 

The numerical results shown in this paper were obtained by 
two alternative methods. Motivated by the success of a sim- 
ilar approach used in random- walks, percolation, and poly- 
mers related problems I23ll24ll25ll26ll . for small systems we 
evaluated the partition function together with the probability 
of every individual state via an exact enumeration method. 
In this approach, first all possible connected configurations 
with a given number of edges are generated successively up 
to M: the graphs with m + 1 edges are constructed from 
the graphs with m edges either by linking a new vertex to 
one of the old vertices, or by linking two previously uncon- 
nected vertices. Next, all possible configurations containing 
M edges are obtained from the combination of smaller con- 
nected graphs, with sizes up to M. Finally, Af a is calculated 
for each topology t a by counting the number of possible per- 
mutations chosen to label the vertices in the state. (For more 
details and a simple example, see appendix A). The probabil- 
ity p a , of a topology t a then can be obtained from 



Pa 



(7) 



Once the set of possible states with appropriate probabilities 
has been constructed, one can evaluate the expectation value 
of any Q thermodynamic quantity using 



(Q) = ]TQ(Up q 



(8) 



The advantage of this technique, beside producing exact 
results, is that the set of J\f a has to be calculated only once, 
independently of the energy functions considered, in contrast 
to Monte-Carlo simulations, where the simulation has to be 
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restarted from the beginning every time we introduce a new 
type of energy. This method is limited by the rapid growth of 
the number of topologies with M. For networks of size seen 
in the real world (M > 10 2 ) the realization of this method is 
clearly unfeasible. 

777. B) Monte-Carlo simulations 

The lattice gas model defined on the edge-dual graph of the 
fully connected network relaxes slowly, because interactions 
are dense Bill and the energy minima are sometimes localized 
in hardly accessible parts of the phase space of the system. A 
good example is the transition from a classical random graph 
- stable at high temperatures - to a star, which is stable at low 
temperatures. The simplest Monte-Carlo rewiring simulation 
(discussed in Sec. II.) tries to move a randomly chosen edge 
to a randomly chosen new location. However, this method is 
very inefficient, if one would like to simulate the condensation 
of edges into a star in a large system. 

There are several simulation tools that can help to achieve 
faster convergence. We have used the so-called parallel tem- 
pering (also called exchange Monte-Carlo) method in several 
cases [27]. Except for first-order transitions, this algorithm 
can be used well to measure the transition at an acceptable 
speed and high precision. 

The algorithm can be viewed as an improved version of 
simulated annealing. Several replicas of the system are simu- 
lated simultaneously, and each of them is connected to a sep- 
arate heat bath. A replica in a hot heat bath will explore the 
"large-scale" structure of phase space, and the motion of a 
replica in a cold heat bath will be restricted to a small part of 
phase space, where it will explore the deep but narrow energy 
wells. 

In the exchange Monte-Carlo method, after a given num- 
ber of conventional update steps within each replica, exchange 
steps are made. Two replicas (with neighboring temperatures) 
are chosen at random, and a Monte-Carlo-type decision is 
made whether the two replicas should be exchanged, i.e., their 
temperatures should be swapped. With Metropolis dynamics, 
if the product of the energy difference between the two repli- 
cas (AE) and the difference of inverse temperatures (A/3) is 
positive, then this exchange is accepted, otherwise it is ac- 
cepted only with probability e A/3AE . That is, a replica with 
a high energy and a low inverse temperature (i.e., high tem- 
perature) will be more likely to remain in its own heat bath, 
whereas a replica with a high energy and a high inverse tem- 
perature (i.e., low temperature) will be likely to be "put" into 
a heat bath with a higher temperature. 

We shall now move on to review some of energy functions, 
that lead to phase transitions, when the temperature is changed 
from zero to infinity at constant (k). Since at T = oo the 
entropically favorable graphs dominate, a minimum require- 
ment for the energy function is that the configurations with 
the lowest energy must also have low entropy. We divide the 
investigated energy functions into two categories. In the first 
category we put the energies that depend on component sizes 



in the graph, the other group contains the single-vertex ener- 
gies. 

IV. CLUSTER ENERGIES 

As mentioned in the introduction, the classical random 
graph model (corresponding to T — > oo) exhibits a phase tran- 
sition when the average degree of the vertices, (k) — 2M/N, 
is varied around (k) = 1. For (k) < 1, the network consists 
of small, disconnected clusters, on the other hand, for (k) > 1 
a giant connected component emerges in the graph collecting 
a finite portion of the edges. Near the critical point the size of 
the giant component scales as ((k) — 1)M. 

Based on the lattice gas analogy we expect that if (k) < 1, 
then for a suitable choice of the energy (one that rewards clus- 
tering) a similar dispersed-compact phase transition occurs at 
a finite temperature T((k)). Such a transition can be best 
monitored by the order parameter <£> s = s max /M (the number 
of edges in the largest connected component s max divided by 
the total number of edges), often used in graph theory |28]. 

The most obvious energy satisfying the above requirement 
is a monotonically decreasing function E = f(s max ). In this 
case the energy is independent of the distribution of the size of 
smaller clusters, or of the structural details of the largest clus- 
ter: only the size of the largest cluster matters. The entropic 
part of the conditional free energy in this case can be estimated 
by counting the number of configurations at given s max . The 
number of different connected configurations of size s max can 
be estimated as s^ax to leading order l29ll . (for an intuitive 
derivation see appendix B). This term has to be multiplied by 
the number of possible selections of these s max vertices out of 
N, which is simply N over s max . The M—s max left-out edges 
can be placed anywhere between the N — s max remaining ver- 
tices, with the restriction that they cannot form clusters larger 
than s max . Since we consider s max = Q S M to be an exten- 
sive quantity and the typical size of the largest component of 
these left-out edges scales slower than M, this constraint can 
be neglected. Hence, the contribution from the left-out edges 



can be well estimated by an (JV — s max ) /2 over (M - 







factor. If we combine these factors together, then in the ther- 
modynamic limit (when TV, M — > oo, (k) = const.) we can 
write 
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$ S MJ\ M-$ a M 



(9) 



Since the energy of the system is a function of the order pa- 
rameter $ s itself, the conditional free energy F(Q S , T) can be 
expressed as 

e -F(* e ,T)/T =A/ - cluge -/(* 8 )/T = e -[/(* s )-Tln^ lu3 ]/^ 10) 

By using Stirling's formula [11 w (l/e) l y/2nl] to approximate 
the factorials and neglecting terms of 0(lnN), for lniV c i us 
we get 
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By replacing 2M/N with (k) in the expression above, the 
resulting conditional free energy can be expressed as 

F($„T) « /($ S M) + MT - 1 - ln«fc»] $ s 

$2 ] 



(Ar)' - 3 (fc) + 2 



(12) 



The simplest choice for an energy function that depends 
only on the size of the largest component is 

/Omax) = "Srnax - -$.M. (13) 

In this case it can be clearly seen from Eq. d!2i that as long as 



T > T c «fc)) 



(fc)-l-ln((fc>)' 



(14) 



the free energy has a minimum at $ s = <i?*(T) = 0, i.e., the 
configuration is dispersed (see main panel of Fig. [2j. When 
the temperature drops below T c ((fc)), the minimum moves 
away from $ s = and a giant component appears. Near the 
critical temperature T c ((fc)), the order parameter at the mini- 
mum of the free energy can be estimated from Eq. d!2t as 



*:(t) 



T^-T-\{k)) 
J (k) 2 -3(k) +2' 



(15) 



indicating that we are dealing with a continuous topological 
phase transition (see inset of Fig. [2}. 

Topological phase transitions of first order are also expected 
to occur for other forms of cluster energies. For example, 
when /(s max ) starts with a zero (or positive) slope. In such a 
case, the behavior of the conditional free energy at very high 
and very low temperatures is similar to the previous case: 
when T — > oo, the energy term /(s max ) can be neglected 
in d!2i . and the entropic term has a minimum at <1> S = 0, 
which is the dispersed state. In contrast, when T — > 0, only 
the energy term remains, resulting a minimum in F($ S ,T) 
at $3 = 1, the compact state. However, there is also an in- 
termediate temperature range, where F($> S ,T) given by (I12> 
starts as an increasing function at $ s = (since the linear 
term in the entropy dominates for small $ s ), then reaches its 
maximum somewhere in the [0, 1] interval and continues as a 
decaying function from that point on (since the higher order 
decaying terms in the energy overcome the increasing terms 
at larger As a consequence, the conditional free energy 
has two competing minima in the [0, 1] interval, a meta-stable 
and a globally stable. The coexistence of stable and meta- 
stable minima at the transition between the phases is a char- 
acteristic of first order phase transitions. A simple function 
of this type is f(s max ) — — s^ ax . For this choice of the en- 
ergy, our numerical results are shown in Fig|3] The hysteresis 




FIG. 2: The phase diagram and the order parameter for the E — 
— Smax energy. Main panel: The white and shaded areas correspond 
to the ordered phase (containing a giant component) and the disor- 
dered phase, respectively, as given by Eq. 1 141 . Inset: The order 
parameter $ = $ s = s max /M obtained from Monte-Carlo simula- 
tions as a function of the inverse temperature for (fc) = 0.1 (trian- 
gles) and (k) =0.5 (circles). Each data point is an ensemble average 
of 10 runs, time averaged between t = 1007V and 5007V Monte- 
Carlo steps. The open and closed symbols represent TV = 500 and 
1,000 vertices, respectively. The critical exponent, in agreement 
with the analytical approximations (solid lines), was found to be 1. 

appearing between cooling and heating supports the theoreti- 
cal considerations about the coexisting minima that indicate a 
first order phase transition, summarized in the inset of Fig|3] 
The analytical conditional free energy gained by substituting 
/(fimax) = — s max mt0 nas a single minimum, when 
T < 80 and when T > 430, in former case at the dispersed 
state (<£> s = 0), in latter case at the connected state ($ s = 1). 
At intermediate temperatures these two minima coexist, pre- 
dicting a first order phase transition somewhere in the mid- 
dle part of this temperature interval. The transition regime 
170 < T < 270 observed in the MC simulation is compat- 
ible with the analytical result, since one does not expect the 
simulation to reveal the meta-stable state beside its dominant 
stable counterpart in case of a very significant difference be- 
tween the depths of the according minima in the free energy. 

We have also investigated the case of /(s max ) = 
—Smax ln(smax), both analytically and numerically. Similarly 
to the previous case, it can be shown that there is a tempera- 
ture regime, where the conditional free energy as a function of 
$ s has two competing minima, hence the transition is of first 
order. 

Since the energy E — f(s max ) depends on a global quan- 
tity (the size of the largest connected component) it might be 
also reasonable to define the energy of the graph as E = 
J2j f( s j)> where the summation goes over each component 
and Sj denotes the number of edges in the jth one. The total 
number of edges, M = J2j s j> i s conserved by the dynamics, 
hence f(sj) must decrease faster than linear to promote com- 
pactification. When a single giant component (containing the 



6 



4> s = Mj / M 

1 i 



i|:lli|:i!,> 



!8 ' • IS 



. 8 a S 



> g 8 » 8 : 



0.5 




0.5 Q 1 



^Il«lfi!!iiflitlf(!ii 



160 180 200 220 240 

T 



* S = M!/M 

l r 

HsSlsSs 



0.5 



•I Mi 



iSilllliillillllliliil 



120 



160 



200 



240 



FIG. 3: If the energy of the graph is E — — s max , then the order 
parameter, $ s = Mi/M, shows a first order transition. (Mi is the 
number of edges in the largest component of the graph.) Each point 
gives the value of $ s averaged between t = 490TV and t = 5007V 
Monte-Carlo steps in a graph started at t = from an Erdos-Renyi 
random graph (o) or a star ( x ). The simulated graph had TV = 500 
vertices and M = 125 edges. The inset shows the behavior of the 
analytical free energy obtained from 1121 using the same parameters. 
The temperature interval in which the two minima (at 3> s = 0, cor- 
responding to the connected state and at <!> s = 1, corresponding to 
the dispersed state) coexist is fully compatible with the numerical 
findings for the transition regime. 

majority of the edges) emerges, its energy /(s max ) dominates 
the energy of the entire graph, and, as a good approximation, 
the above analysis for E = f(s max ) can be repeated, leading 
to first order and continuous phase transitions. For the case of 
E = — Yli=i s i ' our numerical results are presented in Fig|4j 
showing a first order phase transition. 

For the E — — J2iLi s i m ( s i) energy - similarly to the 
case of the E = — s max ln(smax) energy - we found a first- 
order transition between the ordered phase (present at low 
temperatures) and the disordered phase (high temperatures). 
Results are shown in Fig. [5] 



FIG. 4: The order parameter, $ s = Mi/M, for the E = 
— 53 =i s ? energy- (Mi is the number of edges in the largest com- 
ponent and TV C is the number of components in a graph.) Each point 
shows the value of $ s after t = 2007V Monte-Carlo steps in a graph 
started at t = from an Erdos-Renyi random graph (o) or a star 
(x). The simulated graph had TV = 500 vertices and M = 125 
edges. Observe that for intermediate temperatures there are two dis- 
tant groups of states with high stability (at $ s w 0.6 — 0.9 and 
<E> S m — 0.15), and $ s values in the region between these two 
rarely occur, i.e., a first order transition was found. 




V. SINGLE VERTEX ENERGY FUNCTIONS 

Next we turn to another important class of the energy func- 
tions, where the energies are assigned to the vertices rather 
than to the connected components of the graph: 

N 

1=1 

where fcj denotes the degree (number of neighbors) of vertex i. 
This energy is consistent with a dynamics in which the change 
of the degree of a vertex depends only on the structure of the 
graph in its vicinity. The fitness of an individual vertex de- 
pends on its connectivity. The most suitable order parameter 
for this class of graph energy is $ = = k max /M. Again, 



FIG. 5: Distribution of the size of the largest graph component, 
Smax, if the E = — y^._ 1 Si ln(si) energy is used. Main panel. At 
low temperatures the distribution of the largest component has one 
maximum at large values of s max , this is the ordered phase. At high 
temperatures, the largest component is small: the graph is disordered. 
Inset. Using a higher resolution, one can observe that the transition 
from the low-temperature peak (T = 7) to the high-temperature peak 
(T = 9) happens via a bimodal distribution (T = 8, indicated by 
a solid line). This indicates that the conditional free energy of the 
system has two competing minima at the intermediate temperature, 
and the transition is of first order. The graphs used for the simulations 
had TV = 500 vertices and M — 125 edges. Averages were taken for 
10 (main panel) or 200 (inset) simulation runs between simulation 
times of t — 200TV and t = 400TV Monte-Carlo steps using time 
steps of t = TV. 
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due to the conservation of the number of edges, M — J^i ki> 
the single vertex energy f(ki) should decrease faster than 
— ki, if aggregation is to be favored. 

We introduce an alternative form for single vertex energies: 



(17) 



i=l 



where i' runs over all vertices that are neighbors of vertex i. In 
this interpretation, the fitness of an individual vertex depends 
on the connectivities of its neighbors, and vertex i collects an 
energy g(ki') from each of its neighbors. These neighbors 
in turn will all collect g(ki) from vertex i, therefore the total 
contribution to the energy from vertex i is kig(ki). Thus, by 
using 



f(ki) = hg{ki), 



(18) 



the two alternative forms of the single vertex energy, (II 61 and 
(1171 . become equivalent. 

V. A) The energy E — —^kf : mapping to the Ising-model 

A natural choice for the energy of single vertex type is the 
following. Assign the negative energy — J to all pairs of edges 
that share a common vertex at one end. The total energy of a 
given configuration is then 



E 



J 



N 



J 



N 



1 



Y / k l (k l -l) = --J2^ + -JM, (19) 



corresponding to f(ki) = — (J/2)kf, (or equivalently, to 
g(ki) — —(J/2)ki). The constant term in dl9l does not play 
any role in the dynamics, hence it can be omitted. This form 
of the energy is in full analogy with the usual definition of the 
energy 



E = -J 



<a,/3> 



n a n/3 



(20) 



of a lattice gas on the edge-dual graph of the fully connected 
network with nearest neighbor attraction. The summation 
here runs over all adjacent pairs of lattice sites (corresponding 
to possible edges between the vertices of the original graph), 
and n a = 1 if site a is occupied and otherwise. When 
this energy is applied to lattices, we recover the standard lat- 
tice gas model of nucleation of vapors. The negative energy 
unit — J associated with a pair of edges sharing a vertex in the 
original graph is equivalent to the binding energy between the 
corresponding occupied nearest neighbor sites on the edge- 
dual graph. By measuring the energies (and temperature) in 
units of J we can set J = 1, without loosing generality. Thus, 
from now on J will be omitted. 

The lattice gas representation can be further transformed to 
an Ising-model -representation by introducing the s a G [—1,1] 



spin-like variables connected to n a as n a = (1 + s a )/2. The 
energy with the help of the spins is expressed as 

1 iV(JV-l)/2 
E = -J Yj S « S /3"o Yj Sa 



<a,0> 



1 



N(N — 1)(JV- 2), 



(21) 



since the total number of lattice sites equals N(N — l)/2, and 
the number of adjacent pairs of lattice sites is N(N — 1) (N — 
2)/2. This is similar to a ferromagnetic Ising-model in an ex- 
ternal magnetic field. If the number of occupied sites in the 
lattice gas picture equals the number of unoccupied sites, the 
contribution from the external magnetic field vanishes in the 
Ising-model picture. However, in the thermodynamic limit, 
where N — ► oo, M — ► oo, (k) = 2M/N =const., this condi- 
tion cannot be fulfilled, since the total number of sites scales 
as N 2 , whereas the number of particles scales as M with the 
system size. 

For the particular form of f{kj) chosen, the topology with 
the lowest overall energy is a "star" (for simplicity, we con- 
sider M < N), where all the M edges are connected to single 
node. The form of the conditional free energy in this case can 
be estimated as follows. In the $fe > 1/2 regime, where the 
system contains a star of size larger than M/2, the energy of 
this star dominates the rest of the graph. Therefore, in the 
thermodynamic limit, we neglect this latter contribution to the 
total energy and approximate the energy of the graph by that 
of the largest star. Now we estimate the number of possible 
configurations for a given value of In case of a star with 
K = <f>fcM arms, the central vertex can be chosen from TV 
different vertices. Once this is fixed, K edges have to be dis- 
tributed among the N — 1 possible links between the other 
vertices and the central one, yielding a factor of N — 1 over 
K. The rest of the edges that are not part of the star can be 
placed anywhere between the N — 1 (non-central) vertices, 
contributing a factor of (TV - 1)(JV - 2) /2 over (M - K). 

'N — 1\ / (N — 1)(N — 2)/2^ 

K J V M-K 
'N\ f N 2 /2 
K) \M — Kj 



N 



(22) 



Again, we use Stirling's formula to approximate the factorials, 
and neglect terms of the order C(ln N) yielding 



ln7V s t 



M 



N , N 
— In — 

M M 



N 2 N 2 
2M n 2M 



-$ fc ln$ fc -(l-$ fe )ln(l-$ fc ) 

T2 \ f N 2 

$ fc ln 77TT " 1 + $fc 



\ 2M 



(23) 



In the thermodynamic limit, to leading order we receive 

lnJVrtwCaOw-fcfcMlnfJV), (24) 
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1000 

T 



FIG. 6: The order parameter $ = = fc max /M as a function of 
the temperature and the system size for E = ^ . — fef/2 and (fc) = 
0.5. The simulations were started either from a star (corresponding 
to T = 0, solid line) or a classical random graph (T = oo, dashed 
line). Each data point represents a single run, time averaged between 
t = 100iV and 200iV Monte-Carlo steps. The thick solid line shows 
the analytically calculated spinodal Ti = Mj ln(iV). 

where the independent terms were dropped. The resulting 
conditional free energy is expressed as: 

F($ fe ,T) « /(* fc M) + $ fe MT m(AT). (25) 

In the present case, with the /(fej) = — fcf energy, ( I25l > can be 
written as 

F($ fc ,T) « M [-#|M + $ fe Tln(iV)] . (26) 

Note that this approximation would be valid even for < 
1/2, if the energy of the graph was simply defined as E = 

max } ■ 

The parabola given by Eq. \26\ has a maximum at <&k = 
T/M\n{N). When T —> 0, this maximum also shifts to- 
wards zero and T) becomes a descending parabola on 
the [0, 1] interval. This means that the minimum of the free 
energy is at <£>/,. = 1, the star configuration. In contrast, 
when the temperature goes above the Ti = Mj In (TV) spin- 
odal point (thick solid line in Fig.|6j, the maximum moves out 
of the [0, 1] interval and the free energy becomes an ascend- 
ing parabola, resulting in a minimum at a low value of 
(corresponding to an ER random graph). However, this value 
cannot be deduced from Eq. i !26t . because it is a valid approx- 
imation only for > 1/2. For intermediate temperatures the 
maximum of the parabola separates the two extreme topolo- 
gies: the dispersed random graph and the star). One of these 
two extreme states is meta-stable and the other one is stable. 
Due to the limited validity of Eq. ( I26> . the stability of these 
configurations can be studied only for temperatures where the 
maximum of the parabola is well inside the [1 /2, 1] interval. 

The scenario of the transition from a dispersed state to 
the star configuration (see above) indicates that it is a first 
order phase transition. This is well supported by the results 
of both the exact enumeration method and Monte-Carlo 
simulations. For small systems, the conditional free energy 
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FIG. 7: The picture of the conditional free energy at three different 
temperatures for the /(fo) = — fc^ energy, obtained from the exact 
enumeration method plotted together with the prediction of our sim- 
ple theoretical analysis for M = 12, N = 48. At low temperatures 
-F(3 > fe, T) is a descending function on the [0, 1] with a minimum at 
$fc = 1, the star configuration (top figure), on the other hand for 
high temperatures, it becomes ascending for most part, with a min- 
imum at low the dispersed states (bottom picture). There is an 
intermediate temperature regime in between, where the maximum of 
F(&k,T) separates two competing minima (middle figure), hence 
this phase transition is of first order. 

was evaluated via the exact enumeration method for various 
temperatures, and was found to be in qualitative agreement 
with the prediction of the theoretical analysis, as demon- 
strated in Fig. The three different temperature regimes 
described in the previous paragraph can be recognized in the 
behavior of the exact T($fc, T) as well. Furthermore, in the 
intermediate temperature regime, where the conditional free 
energy has two competing minima, the spinodal curve can 
also be constructed as is shown in Fig. [S] For large enough 
systems, in MC simulations a sudden change of the order 
parameter between zero and one can be observed as shown in 
Fig lU The hysteresis appearing between cooling and heating 
is consistent with a first order transition. 

V. B) The energy E = — ^fcilnfci : continuous phase- 
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FIG. 8: The spinodal curve obtained from the exact enumeration 
method with E = - ^ kf, for M = 12, TV = 48. At low and 
high temperatures, the conditional free energy F(&k, T) has a sin- 
gle minimum (plotted with squares). At intermediate temperatures 
(in between the two dotted lines) there are two competing minima. 
In this latter temperature regime, the spinodal curve is obtained by 
plotting the maximum of F(<f>k,T) (represented by stars), besides 
the two minima. 

transition 

Another application-motivated choice for the single ver- 
tex energy is /(fcj) = — fcjln(fcj), or equivalently, g(ki) — 
— ln(fcj), inspired, in part, by the logarithmic law of sensa- 
tion. It is the logarithm of the degree of a vertex that its 
neighbors can sense and benefit from. In this case the con- 
figuration with the lowest energy is a fully connected sub- 
graph [or almost fully connected if M cannot be expressed 
as n(n — l)/2]. On the other hand, the star configuration 
is also quite favorable, since the energy of both the maximal 
possible star and of the maximal possible fully connected sub- 
graph scales as — M In M to leading order. Amongst the sub- 
dominant terms in the energy, there is a difference in the order 
of ^/Mln^/M between the two, in favor of the fully con- 
nected subgraph. As before, we choose the order parameter 
to be $ = $/j = fc max /M, since this can easily distinguish 
between these two configurations: fc max y/2M for a fully 
connected subgraph counting M edges, and fc max w M for a 
star. 

Our MC simulations demonstrate (Fig. |9} that as we cool 
down the system, first the edges of the dispersed random graph 
assemble to form a configuration with a few large stars (shar- 
ing most of their neighbors), and then at lower temperatures 
the graph is rearranged into an almost fully connected sub- 
graph. This is consistent with the fact that beside the slight 
energetical disadvantage, the star configuration is entropically 
more favorable when compared to the fully connected sub- 
graph; therefore the latter configuration can take over only at 
very low temperatures. The hysteresis near the few large star 
vs. fully connected subgraph transition suggests that it is a 
first order phase transition. On the other hand, the transition 
between the dispersed state and the few large stars is accom- 
panied by a singularity in the heat capacity (also seen with 
the exact enumeration method), and no hysteresis is observed, 




FIG. 9: Phases of the graph when the energy is E = 
— ki ln(fc). (a) The largest degree fc max for N = 10, 224 ver- 
tices and M = 2, 556 edges. Each data point represents a single run, 
time averaged between t = 5, 000/V and 20, 000JV MC steps. The 
data points are connected to guide the eye. There is a sharp, continu- 
ous transition near T — 0.85 and a first-order transition (with a hys- 
teresis) around T — 0.5 — 0.6. (b) The three different plateaus in (a) 
correspond to distinct topological phases: fc max = 0(1) to the clas- 
sical random graph, fc max = O(M) to the star phase (a small number 
of stars sharing most of their neighbors) and fcmax = 0(vM) to the 
fully connected subgraph, (c) The (cumulative) degree distribution 
at T = 0.84 and t = 600/V follows a power law. This shows that the 
degree distribution decays as a power-law with the exponent 7 « 3. 

indicating that it is a continuous phase transition. 

For $fc > 1/2 Eq. d25i can be used again as a good approx- 
imation for the free energy of the graph, since the compact 
cluster arising from the dispersed state is rather star like. By 
plugging f($ k M) = -($ fc M) ln($ k M) into that expres- 
sion, we get 



F($fc,T)wM(r-l)ln(JV)#fc 



(27) 



to leading order, which is linear in In agreement with our 
observations above, this formula predicts that for T < 1 the 
star is a stable configuration (<£>A;=1 is a minimum of the free 
energy), and for T > 1 it becomes unstable. The transition at 
T = T c = 1 is thus step-like with no hysteresis, indicating 
a continuous phase transition with an infinitely large critical 
exponent. We assume that the observed deviation of T c from 
1 in the MC simulations is a finite size effect. 

V.C) Relation to growth with preferential attachment 

A remarkable feature of the MC dynamics is that in case of 
the energy f(ki) = —kiln hi, by crossing T c from above, 
a scale-free graph (with a degree distribution ~ k 1 with 
7 ~ 3) appears at some point of the evolution of the graph 
from the random configuration towards the star. This supports 
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the notion that scale-free graphs are temporary (dynamical) 
configurations, not typical in equilibrium distributions. The 
MC dynamics is governed by the change of the energy asso- 
ciated with the reallocation of an edge. Estimating the energy 
change of a vertex by the derivative of the single vertex energy 
f(ki) = ~ki ln(fci), we get AE = 1 - ln(fci). Plugging this 
into the Boltzmann factor, exp[-AE/T], at T = T c = 1 we 
get a quantity proportional to fej for the acceptation/rejection 
ratio of a randomly selected move. Since the preferential at- 
tachment in the Barabasi- Albert model |4] is proportional to 
fej, it is natural that our dynamics also produces scale-free 
graphs. 

Another interesting aspect of the /(fcj) = —kiln fc, energy 
is that the configurations in the two compact phases resemble 
the two major graph topologies obtained in Ref. frsll . by 
optimizing the network for local search with congestion. Our 
intermediate phase with a few large central hubs sharing 
neighbors is similar to the optimal topology for a small 
number of parallel searches, whereas the low-temperature 
configuration, the fully connected subgraph resembles the 
homogeneous topology optimal for a large number of parallel 
searches. However, an important difference between the 
two problems is that in our case a vertex is allowed to lose 
all of its connections under the restructuring process. The 
two "similar-to-optimum" configurations appear as a natural 
consequence of the underlying dynamics. This observation 
suggests a potential application of the presented theory: 
tackling problems related to graph topology optimization by 
simulated annealing techniques. 



VI. THE GRAND CANONICAL ENSEMBLE 

In the grand canonical ensemble, the degree distribution can 
be expressed as frill 

P k = C . (28) 

where C is a normalization factor and the chemical potential 
\i is adjusted to give the correct (fc). For f(k) — — fcln(fc), 
using Stirling's formula, the distribution takes the form 

e -0-i)fc 

P k = C — ^^k {1/T - 1)k . (29) 
V27rfc 

When T > 1, this has a tail, which decays faster than ex- 
ponential, consequently, each vertex has a small degree. For 
T < 1, on the other hand, the tail becomes divergent, signal- 
ing a phase transition at T = T c = 1. Note however that in 
the T < 1 temperature range, due to the non-extensive contri- 
bution of the diverging degrees, the ensembles are not equiv- 
alent, and the grand canonical description loses its validity. 

At the critical temperature, the grand canonical description 
might still be valid. Choosing a more general single ver- 
tex energy, f(ki) = —(hi — a)\n(ki), and setting (k) such 
that \x — 1, the degree distribution acquires a power law tail 
(Pk ~ fc~' Q+1 ^ 2 - ) ) and the network becomes scale-free at this 
temperature. We have to stress though that the scale-free net- 
work at T c is not general: for // > 1 the tail decays exponen- 
tially, and for /i < 1 the tail diverges. 

VII. SUMMARY 



YD) Topology-dependent non-extensiveness of the energy 



Both types of the single vertex energy functions discussed 
in the present section lead to compact configurations at low 
temperatures, for which the most highly connected vertices 
possess macroscopic numbers of edges. As a consequence, 
the energy of the system scales differently with system size 
at high and low temperatures, and diverges differently as 
N — > oo. At high temperature, the system consists of many 
small unlinked clusters of about the same size, therefore a 
change in the total system size affects only the number of the 
clusters, and the energy scales as N. On the other hand, when 
f(ki ) = —ki ln(fcj), at low temperatures the energy of the star 
and the fully connected subgraph scales as Nln(N); in case 
of f(ki) = — jfc?, the energy of the star scales as N 2 . Thus 
(unlike, i.e., in the mean-field Ising model), there is no way 
to choose an appropriate coupling constant that could render 
the energy extensive in all topological states simultaneously. 

Nevertheless, the dispersed state (having an extensive graph 
energy) can equally be studied in the grand canonical ensem- 
ble. 



We studied the restructuring in networks using a canoni- 
cal ensemble, where temperature corresponds to the level of 
noise in real systems and the energy associated with the dif- 
ferent configurations accounts for the advantage gained or lost 
during the rewiring of the edges. We found that for vari- 
ous types of energies, first order and continuous phase tran- 
sitions may appear when changing temperatures. In case of 
the E = — s ma x energy, if (k) < 1, a dispersed-loose phase 
transition occurs at a finite temperature, equivalent to the per- 
colation phase transition of classical random graphs when (k) 
is varied around (k) = 1. We obtained a simple expression 
for the T c ((fc)) critical line separating the two phases in the 
[(k) , T] plane from a theoretical analysis of the conditional 
free energy. For other forms of the energy depending on the 
size of the largest cluster we found first order phase transi- 
tions. We also studied the effects of different single vertex 
energies, namely the E = — kf and E = — J2i &i Wh) 
cases. The network in the former case exhibits a first or- 
der phase transition from a dispersed state to a star-like state, 
where nearly all edges are linked to a single vertex. With the 
— J2i ki ln(fci) energy, the dispersed state transforms into a 
compact one with a few large stars via a continuous phase 
transition. In the critical point, scale-free networks can be re- 
covered. At lower temperatures another transition occurs (this 
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time of first order), where the configuration is turned into a 
fully connected subgraph. 

Although in this paper we assumed that (k) < 1, this is not 
a necessary requirement, when the energy is assigned to indi- 
vidual vertices. For large average degree ((fc) > 2) the only 
difference is that one vertex cannot collect all the edges, and 
thus, several stars appear in the "star" configuration. Further 
interesting directions in the context of the above study include 
the investigation of additional relevant forms for the energy 
[e.g., E — (k — n) 2 with n > 1.5] and the joint effects of 
restructuring and growth. 
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APPENDIX A 

In the exact enumeration method, as mentioned in sec. III. A, 
the first step is to generate all connected graphs with m + 1 
edges from the connected graphs with m edges, either by 
connecting a new vertex to the core or by introducing a new 
link. In order to avoid double counting, every new graph ob- 
tained this way is compared one by one to all already revealed 
topologies using the following algorithm. Two graphs of iden- 
tical topology have identical degree distribution also, there- 
fore this property is checked first. In case of perfect match, 
the vertices in both graphs are labeled in such a way, that a 
given index belongs to vertices with equal number of links in 
the two graphs. Next, for each index in one graph, the set of 
the neighbors indices is compared to its equivalent index set in 
the other graph. If not all sets are identical, then the labels in 
one of the graphs have to be permuted until perfect match be- 
tween the neighboring relations is reached. (Obviously, labels 
are interchanged between vertices of same degree only). If the 
perfect match in the neighboring relations cannot be achieved 
for any permutation of the indices, the two graphs are of dif- 
ferent topology. 

When a new topology is obtained, the corresponding com- 
binatorial factors can be generated in a similar manner, by 
counting the number of permutations of the indices in the 
graph that lead to the same neighboring relations (same neigh- 
boring index sets) as the original indexing. 

As a simple example, we demonstrate the evaluation of J\f a 
for all states in case of M — 3, TV > 6. The construction of 
the connected graphs, and the possible topologies are shown 
below: 
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1) 2) 




replaced by the following graphs : 



In case of a topology that does not possess any symmetries, 
JV is simply A!/(A — N t )l, where N t is the number of ver- 
tices included in the topology. In general this initial Af has 
to be further divided by the number of those permutations of 
the indices of the vertices that leave the topology unchanged. 
Therefore, if the topology contains n identical subgraphs (like 
in case of state a = 1 above, where the topology is built up 
from three identical subgraphs) the initial value of J\f has to be 
divided by n\. Furthermore, if any subgraph in the topology 
remains unchanged for I permutations of the indices within it- 
self, J\f has to be divided by I. In the example above, for the 
states a = 1 and a — 2, for all subgraphs, I = 2, in case of 
the states a = 3 and a = 5, I = 3!, and for the state a = 4, 
I = 2. 

Altogether, in the chosen case, the J\f of the five possible 
states can be expressed as 



AH 



A! 



(A-6)!2 3 3! 
A' 

M= (A34)!3l> * = 



(A-5)!2 2 ' 

■Ms = 



A! 



(A-4)!2' 



(A — 3)!3! ' 



To provide a simple example of an application, we show 
the first few most probable states in case of E = — ^ ki In ki 
at T = 0.65 : 



[i=0 W)2 1 d p=(MXV)l(j4 p=U0UM46 p=()(H}6l77 p=0.0O5798 p=0.005495 p=(MM)5.«.? p=0 i)I)5:dp 



p=O.005O97 p=0.005026 p=0.004NHO p=() (Ki4(.l)^ p=(i |)((4557 p=() O044P7 p=() (I(t4.i70 p=0.004321 



When the temperature is lowered to T = 0.3, these are 



APPENDIX B 

For simplicity, we shall consider tree like clusters only and 
neglect the clusters with loops. Since the chances of a com- 
ponent containing a closed loop of edges goes as A -1 when 
(k) < 1 and no giant connected component can be found in 
the system, this is a valid approximation in the thermody- 
namic limit The number of possible trees of size s in 
an undirected network can be estimated as follows. We pick 
a random realization of a tree sized s (meaning s edges and 
s + 1 vertices), and we choose a vertex in it to be the "root" of 
the tree. Starting from this root, we descend through all pos- 
sible paths until we reach all the branches, and on the way we 
replace the undirected edges with directed ones pointing from 
the vertex closer to the root towards the vertex farther away 
from the root. This procedure results in a directed tree, where 
each vertex (except the root) has one and only one incoming 
edge and n > outgoing edges. Then, another realization 
of a tree can be obtained from the present one by choosing a 
vertex, and moving the other end of the incoming edge from 
its original place to a new vertex. Of course, this new ver- 
tex cannot be one of the "descendants" of the selected vertex, 
since that way we would create a loop and split the tree into 
two unconnected parts. Nevertheless, if s is large enough, for 
the majority of the vertices this restriction eliminates only a 
negligible part of the possible rewirings. Therefore we may 
estimate the number of possible new trees obtained from the 
rewiring of the incoming edge of a single vertex by s, and the 
total number of trees of size s by s s . 



